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I. INTRODUCTION 



The Lattice Boltzmann method (LBM) is a useful tool for simulations of the complex 
hydrodynamic phenomena such as turbulent flows, multiphase flows and suspensions. It 
is even believed that the subsequent development of the lattice Boltzmann method may 
provide a new paradigm in the kinetic modeling due to its mathematical simplicity and 
computational efficiency 

One important issue, which attracted considerable attention recently, is the enhancement 
of stability of the method ^. We remind that the early predecessor of the LBM, the lattice 
gas method of Frisch, Hasslacher and Pomeau pj, was consistent with the H theorem by the 
microscopic detailed balance and was supported by nonpolynomial equilibria (maximizers 
of the Fermi-Dirac entropy). The method was was unconditionally stable, but in the course 
of the subsequent development of the LBM this feature of the unconditional stability, which 
otherwise would distinguish the lattice Boltzmann method among other methods of compu- 
tational fluid dynamics, was gradually lost. The reason why this happened can be traced 
back to the earliest versions of the LBM, derived from the lattice gas model, where the true 
nonpolynomial equilibria were replaced by their low-order polynomial approximations. Of 
course, this was motivated by a search for computationally more effective schemes, begin- 
ning with the work of Higuera, Succi and Benzi M, and which eventually culminated in the 

fin 

athermal single relaxation time lattice Bhatnagar-Gross-Krook model (LBGK) p, |6|] based 
on polynomial equilibria. Later studies aimed at restoring the if -theorem for the LBM are 
well documented (see, e. g. j^J). At present, the entropic LBGK models which combine 
computational efficiency of the standard LBGK with the unconditional stability pertinent 
to genuine kinetic models are constructed [7[ . First results of simulation of high Reynolds 
number flows confirm the theoretically expected significant overall gain in performance of 



the LBM by using the entropic formulations 

The present study is motivated in part by a recent publication entitled "Nonexis- 
tence of H theorems for the athermal lattice Boltzmann models with polynomial equilibria" . 
Therein, the authors demonstrated that polynomial equilibria used in the lattice Boltzmann 
method on the two-dimensional triangular lattice (the D2Q7 model, see section IIII|) do not 
minimize any convex function, at least when parameters of these local equilibria are in a cer- 
tain range. Since the fact that polynomial equilibria are at odds with the maximum entropy 
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principle (and thus they are not local equilibria in a thermodynamic sense) was pointed 
out already for some time (see, for instance, lp|), and actually the experience gained from 



many studies of various discrete velocity models does not indicate that polynomials are 
expected as local equilibria, it is therefore not surprising that the computation jsj failed to 
derive a suitable entropy function for the D2Q7 model. 

In this paper, we revisit the problem of finding the H function for the D2Q7 model. 
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Our approach to derivation of the H function follows the method suggested earlier in 
The straightforward computation presented in detail results in the unique Boltzmann-like 
entropy function for the D2Q7 setup, thereby enabling construction of the entropic lattice 
Boltzmann models for this lattice. 

The structure of the paper is as follows: In section |H1 for the sake of completeness, we re- 
view briefly the athermal LBGK and the athermal entropic LBM, in particular, the entropic 
LBGK. This section contains no new results. In section IIII[ we consider the 2DQ7 lattice, 
and find the appropriate H function. In section \W\ we find an approximate solution to the 
equilibrium distribution and discuss its accuracy. Finally, in section El we put the entropic 
lattice Boltzmann method into a perspective with other recent approaches to stabilization 
of the athermal lattice Boltzmann models. 



II. LATTICE BOLTZMANN AND ENTROPIC LATTICE BOLTZMANN 

In this section, for the sake of completeness, we briefly compare the entropic lattice 
Boltzmann method (ELBM) and the standard lattice Boltzmann method (LBM) of athermal 
hydrodynamics. 

In both the LBM and ELBM methods, one considers populations /j of discrete velocities 
Cj, where i — 1, . . . , m, at discrete time t. The discrete velocities form the links of a regular 
and sufficiently isotropic lattice, and it may also include a zero vector. It is convenient to 
introduce m-dimensional population vectors /. In the isothermal case, local hydrodynamic 
variables (density p and momentum density pu) are defined at lattice sites r as: 

in 

1=1 (i) 

pu = ^2dfi(r,t). 
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The ELBM begins with finding a convex function of populations, H, which satisfies the 
following condition: If f^ q (p,pu) minimizes H subject to the hydrodynamic constraints (JTJ), 
then / eq also verifies the Galilean invariance of the stress tensor: 

m 

^ C ia C if3 f- q (p, pu) = p c 2 s 5 a)3 + pu a up. (2) 

i=l 

Here c s is sound speed. 

The H function which satisfies this condition to the accuracy of u A , and thus is valid 

to all purposes of incompressible simulations was derived in [12J for the the D1Q3 and the 
D2Q9 lattices. Later, this result was extended to the three-dimensional D3Q27 lattice in 
I3I ] . Recently, other H functions which verify equation (J2J) to the accuracy of w 4 were found 

for isotropic Bravais lattices in [Tj|, in particular for the D2Q6 model. In order to illustrate 



;his, we list 



7. 
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rere the results for the H functions and their minimizers for the DkQ3 k lattices 



3, 15, 16,3. be the spatial dimension. For D = 1, the three discrete velocities 

are C = { — 1,0, 1}. In higher dimensions, the discrete velocities are tensor products of the 
discrete velocities of these one-dimensional velocities. The H function is Boltzmann-like: 

* = £>ln(A). (3) 

1=1 N ' 

Here w\ is the weight associated with the ith discrete velocity Cj. For D — 1, the weights 
corresponding to the velocities C = { — 1,0, 1} are w = {|, |, |}. For D > 1, the weights 
are constructed by multiplying the weights associated with each component direction. 

The local equilibrium which minimizes subject to the fixed density and momentum 
reads: 



a=l 

The speed of sound, c s , in this model is l/\/3. 

Once the entropy function H is found, the basic equation of ELBM to be constructed 
and to be solved is 

fi{r + d, t + 1) - f { (v, t) = -Pa[f(v, t)]Ai[f{r, t)\, (5) 

where the right hand side represents the collision process. The m-dimensional vector function 
A (so-called bare collision integral), must satisfy the conditions: 



^jAj{l,Cj} = (local conservation laws), 



i=l 



m dH 

a = < (entropy production inequality). 

Moreover, the local equilibrium vector / eq must be the only zero point of A, that is, 
A(/ cq ) = 0, and, finally, / eq must be the only zero point of the local entropy produc- 
tion, <x(/ cq ) = 0. 

The conditions just listed are standard requirements taken directly from the well known 
theory of the continuous Boltzmann equation. The lattice specifics comes through the factor 
f3a[f(r,t)] in equation (JSJ). Here (3 is a fixed parameter in the interval [0, 1], and is related 
to the viscosity (the limit j3 — > 1 is the zero viscosity limit). The scalar function a is the 



nontrivial root of the nonlinear equation |12L Il4l Il8| | : 



H(f) = H(f + aA[f]). (6) 

It is the function a which ensures the discrete-time i?-theorem, unlike in the continuous-time 
case where it essentially suffices to ensure only the entropy production inequality. Function 
a has to be computed numerically on each lattice site at each time step. In the fully resolved 
hydrodynamic limit, when / — > f eq , the solution a(f) tends to its limiting value a = 2. 

In practice, construction of the bare collision integral A is guided by simplicity. For the 
H functions (JHJ), the local equilibria are given by explicit formula (jlj), and thus the BGK 
form, A = / — f eq becomes available for efficient numerical realizations. We further refer to 
this model as the entropic lattice BGK (ELBGK). A gradient single- relaxation time models 
circumventing the BGK form, and which are readily constructed once just the if-function 
is known, were developed in |l3L Il5l Il9| (see also their discussion in the context of reaction 
kinetics, \2^). 

The standard (nonentropic, second-order polynomial) LBM can be considered as a trun- 
cation of the ELBM's just discussed. This truncation is done in three steps. First, the local 
equilibria are replaced by their second-order in u polynomials. For example, for D = 2, if 
we expand the local equilibrium (JU) to second order in u/c s , we derive the polynomial equi- 
librium of the standard D2Q9 model In order to distinguish between the local equilibria 
when they are minimizers of appropriate entropy functions and the kth order polynomial 
approximations to them, we denote the latter as / . Second, instead of nonlinear bare 
collision integrals one considers linearized forms, Aj = Y^jLi^ijifj — / )■ ^ ne simplest 

~(2) 

option is the BGK form which becomes always available, Aj = /,■ — f) . Third, when the 
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latter expression is substituted into the right hand side of equation (jSJ), the root of the 
equation © is replaced by a = 2. Obviously, these operations do not leave the if-theorem 
intact, and it is maybe not surprising that a polynomial cannot be itself a minimizer of any 
entropy function anymore, as was demonstrated for the D2Q7 lattice in j^. In the subse- 
quent section we find the H function for this lattice without assuming a polynomial ansatz 
for the equilibrium. 



III. H FUNCTION FOR THE 2DQ7 MODEL 

The discrete velocities of the D2Q7 model at each site of a planar triangular lattice 
consist of a zero vector Co = 0, and of six vectors of equal length, Cj, where % = 1, . . . , 6, 
Cj = (cos [ii — l)7r/3), sin ((i — l)7r/3)). The explicit form of the seven-dimensional vectors 
corresponding to the x and y components of velocities are as follows: 



C x = {0, 1, 1/2, -1/2, -1, -1/2, 1/2} , (7) 
/3 

C y = ^j- {0,0, 1,1, 0,-1,-1}. (8) 
Accordingly, the population vector / is, 

/ — {/(b /i) /2j /3j fi-, hi fe}- (9) 
By symmetry arguments, it is sufficient to seek the H function of the form, 

6 

H(f) = h (f )+J2Hfi), (10) 

i=l 

where ho(x) and h(x) are two convex functions of one variable to be determined. [We could 
equally begin with a more general form, H = 5^i=o^*(/i)' assuming a separate unknown 
function for each population. However, the result would be the same.] 

The local equilibrium / eq is the minimizer of the function H (fTU|) subject to the con- 
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stramts, 



ft > 

E /r<^ = /r - /r + 5 (/ 2 eq + / 6 eq - - / 5 eq ) , (ii) 



PUx = 

i=0 

6 nz 



i=0 

Let us denote the inverse of the derivative of h and /i as = [^<o]~ anc l A 4 = WY~ ■ 
Then, the formal solution to the minimization problem reads, 

fo* = Po (x) , 

/r = /i(x+o, 

fit* = P \X + ^Cx + Cy 

fP =fi(x~ \Cx + Cy J , (12) 



/T = vix- Cx) 



fP =fi\X - \Cx - Cy I , 
fp = P (x + \Cx - Cy 

Here Xi Cx and Cy are Lagrange multipliers associated with the constraints. 

By choosing various pairs of functions for /x and /i constitutive relations for the stress 
tensors (J2J) is satisfied with varying degrees of accuracy, and the goal is to find such functions 
/io and \i for which the error reduces to higher orders in the powers of velocity u. We rewrite 
the constitutive relation (j2J) in the form of a discrepancy of the components of the stress 
tensor as: 

r - nr 2 + ( pUx } 2 _ ma r 

J- xx — P^ s i / jJi ^tx^ixi 

1 i=0 

Tyy = PC 2 S + — ^ E fi^CiyCiy, (13) 

P i=0 

_ {pu x ){pu y ) _ f eq r r 

xy s_j Ji ^ix^iy 

P i=0 
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Note that the sound speed is not yet defined in these expressions. In fact, as we will see 
soon, the choice of the sound speed is a solvability condition in our procedure. 

The next step is crucial |l2j: We are going to find such functions /j,q and [i which contain 
no discrepancy up to the orders Q, and ( x ( y . In order to do this, we expand the terms 
in equations (f*T3|) to relevant orders around the point ( x = ( y = 0: 

p(x,C) = Po(x) + QKx) + l^(x)£ + ^'(x)C 2 y + 0(\t\), 
pu x ( X ,C) = 3^(x)C, + 0(|C| 2 ), 

pu y {xX) = 2v / 3 / u'(x)C, + 0(|C|), 

i2fr(xX)c ix c ix = Mx) + -/ix) (^Ic + K) +o(\C\ 2 ), 
j^n\xX)c iy c iy = Mx) + ~/{x) Qc 2 + 2C 2 ) +o(KI 2 ), 
£ r (x, C)c ix c iy = ^"{ x )Uy + o(|CI 2 ). (14) 

Here primes denote corresponding derivatives. 

Substituting expansions (|T4|) in equations (fT3j) . we require that discrepancy (fT3j) vanishes 
to second order in 

At zero order we have two identical equations: 

T<g = c 2 (/io(x) + 6/x(x)) - MX) = 0, 

= c 2 (/i (x) + 6/i( X )) - 3KX) = 0, (15) 

whereupon, 



C 2 = W . (16) 

Mx) + 6 Mx) 

There are no terms linear in £ in the expansion of T a p. At second order we get the following 
equations: 

is - { (§<* - 1) + c + - 1) / we, in) 
*s> - (I^ 2 - 1) + { (2^ - §) „»w + ^^fe} c;, us) 
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3 c s 2 -|U"(x)=0, (22) 



(2) 

We now require T v i — independently of the values of £. Thus, setting to zero each term 
in front of each combination CaCp i n equations (JT7J), (fTH|) . and (fT^j). we obtain five equations: 

2c 2 - i) = 0, (21) 

r 2 - 

)™*-,"(x)=0- (24) 

mx) + 6//(x) 

Equations (J2~T|) and (j22j) ar e identical. Assuming //' 7^ 0, equations (j2*Tj) and (j22j) fix the 
value of sound speed: 

c s = \ (25) 

It is straightforward to demonstrate that with the value of sound speed (|23j) the remaining 
three equations are resolvable. Indeed, substituting (|25)l into the zero-order relation (|16jl. 
we obtain, 

M = 6/i(x)- (26) 
Substituting equations (|23|) and (|2l)|) into equations (|2H|). fl2HJ), and (j2U), we find out that 
each of the latter three equations reduce to the same ordinary differential equation: 

iv'(x)} 2 



Mx) 

The general solution to this equation is 



/Ax) = 0. (27) 



fjL(x) = Aexp(x)+B. (28) 

Now, from the requirement that fi" ^ we get A 7^ 0, and furthermore A > by required 
concavity of H. Substituting the general solution (J2S|) into the equation ()27|). we find B = 0. 
From equation (|2l)j) we then have fio(x) = 6Aexp (%). Therefore, h'{x) = /i _1 ( x ) = In (x/A), 
so that /i(x) = x(hx(x/A) — 1) + fci, and similarly, /io(^) — x(ln(x/(Q A)) — 1) + &2, where k\ 
and /C2 are arbitrary constants. Thus, we obtain the family of Boltzmann-like H function of 
the 2DQ7 model: 

^KK^H+pMCH^ (29) 
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As is well known, adding a linear combination of the locally conserved quantities to the 
entropy function is immaterial, so we can fix A = 1/e, where e is the base of natural 
logarithm, and C = 0: 

H = fo\n(^+Y^fi^(fi)- (30) 
^ ' i=i 

In the next section, we shall study the equilibria corresponding to the entropy function ([3~U|).. 



IV. EQUILIBRIUM POPULATIONS 



By construction, the expansion of the minimizers of H (J3*0"j) around the point u = to 
the order u 2 satisfy all the usual requirements needed to derive the athermal Navier-Stokes 
equations (cf. Ref. 12jJ, see also below for the present case). The exact minimizers are 
nonpolynomial, and are not always available in a closed form. Nevertheless, a glimpse of the 
full solution is possible since the explicit solution can be computed for a few special cases. 
For example, when u x = \^3u y , we find the exact minimizer of the H function (}3*0"j) . 



eq 




P 



1 



1^(9 + 48^.) 
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fT 

req 



r eq _ r 



P_ 

36 

req 
feq JO 



2^/(9 + 48«2) _ 3 _ \2u x 



(31) 



eq 



6 

2pu x 



+ fT- 



The exact solution (}3Tj) will be used below to test the accuracy of various approximations 
to the equilibrium populations. 

For practical realizations, we describe a systematic procedure to obtain the equilibrium 
in a series representation. The procedure relies on the fact that at zero velocity equilibrium 
is known exactly, 



f?(p, 0) = /f } = pwi, w = 1/2, Wj = 1/12, j = 1, 



6. 



(32) 



Once the exact solution for zero velocity is known, extension to u ^ is found by pertur- 
bation. Specifically, the Lagrange multipliers are expanded as 



X 



(33) 



n=0 



n=0 
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where we have introduced a bookkeeping parameter e, , such that u a — > eu a , and e is set 
to one in the end of computation. This series representation of Lagrange multipliers, when 
substituted into equations for the constraints, induces polynomial approximations f- k ^ of 
increasingly higher order, 

k 

/7*> = 5y^>. (34) 

n=0 

Further, we seek the expansion parameters consistent with the conservation constraints at 
all orders. This translates into a set of linear equation solved recursively: 



Eif^' Ecwf^p^, k>2. 



i=0 



For example, at first order of this expansion 

E /i 0) + e ] a x + c^c iy \ = o, 
i=i 

E/i 0) ^ [x (1) + c^ + c^] = p^, 
i=i 

Eif^ [x (1) + e ) ^ + c«a,] = p« r 



i=i 



By solving this linear system of three equations we get, 



= 0, 



Similary, at second order: 



(0) 



1=1 



i=l 



2 X w + 2(d 2) a a ) + (d 1) c ta y 

2 
2 

2 X W + 2(d 2) C7 fa ) + (d 1) C7 ia ) : 



= 0, 



= o, 



= 0, 



(35) 



(36) 



(37) 



(38) 



which gives: 
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X ^ = -2(ul + ul), Ci 2} = 0, Cf = 0. (39) 

Solutions at higher orders are easily found using symbolic computation facilities. In Table 
HJ we present terms of the expansion of the Lagrange multipliers f!33|) up to 8th order in u. 



TABLE I: Expansion coefficients of Lagrange Multipliers. 



k 




<^x 


Jk) 


1 





Au x 


4 % 


2 


-2 {ul + ul) 








CO | 











4 











5 





^X + ^ 


^Uy (5^+3^) 


6 


-§(2*4 + 15*4 u 2 y + Su 6 y ) 
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-§u y (Uul + 21uiu 2 v + 9u 6 y ) 


-§(5*4 + 42*4*4 + 21*4 
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§(5*4 + 56*4 u 2 y + 42u x *4 + 9t^) 









This simple procedure gives the equilibrium distribution at eventually any desired order 
of accuracy. For example, the 0(u 6 ) accurate approximation of the equilibrium is, 




We note in passing that, by construction, each approximate equilibrium population, 
\p, u), satisfies exactly the consistency relations, p{~f^ \p, u)) = p, and pu(f^ \p, u)) = 
pu, at each order k. 

With the above results, we proceed now to evaluate the higher order moments of the local 
equilibrium pertinent to establishing the hydrodynamic limit of the model. The components 
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of the equilibrium pressure tensor are: 

p:i = p (c 2 s + Ux 2 ) - P 



u 



2 2 (u x 2 + 3u y 2 ) 



+0(u b 



Pp y = p(c 2 + u 2 ) - p 



(u x 4 + 3 uJ) 



+0(u e 



(41) 



4 

= PU X Uy ~ -pU^Uy+OiU 6 ), 

4 

Pyx = pU y U x - ~pU x 3 U y +0(U 6 ). 

Here the leading in u terms are responsible for Galilean invariance of the momentum equa- 
tion, whereas the under-braced term give the leading-order deviations. All these deviations 
are of order 0(u A ), as expected by construction of the H function (|3T)j) . We note in passing 
that in the D2Q7 case under consideration, the accuracy of the equilibrium stress tensor, 
as compared with the entropic 2DQ9 model (j3J), is reduced on two counts: First, since c 2 is 
smaller, the effect of deviations is larger, and secondly, the off-diagonal part in equation ()41)1 
is 0(u 4 ) accurate whereas it is 0(w 6 ) accurate in the entropic 2DQ9 model (@J). Similarly, 
the contracted third order moment, Q^ga = Y^=o c iaC 2 fi q , related to the equilibrium energy 



flux, is 



Q7w = (p c s( d + 2 ) + p^ 2 K - pu a u 2 



(42) 



Again, the first term in the latter equation correspond to the well-known result of the 
continuous kinetic theory, whereas the under-braced term is the deviation. Due to the 
lattice symmetry, the moment Q^bb ls ^ ne same f° r an Y equilibrium on the D2Q7 lattice, 
and it is only accurate to the linear order. As we will see it below, it is the accuracy of Q eq 33 
which dictates the working window of the method. 

With the expansion method described above, one can develop the D2Q7 ELBGK model. 
The bare collision integral (cf. sectioning is assumed in the form, A, = — (/» — ), where it 
should be decided that uo to which order k is appropriate. Three condition guide the choice 
of the order A;: 

• Approximation of the populations should be good enough to enable solving for the 
entropy estimate 
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• Deviations in the stress tensor and in the energy flux should be small. Since the 
deviations in the stress tensor are one order higher as compared to the energy flow, 
the latter will be most crucial; 

• The velocity window should not be too close to zero in order to avoid large computa- 
tional time, and also in order to avoid computations with too small numbers. 

In order to establish the working window and the required order of accuracy of the 

~ CO 

approximation / , we test the results of the expansion against the exact solution (|31|). In 
Fig. Q zero velocity component /q is compared with the exact solution for various k. It 

~(2) 

is clear that the usual second-order approximation /q is good enough only for u < 0.001. 
However, for a larger velocity window, u < 0.1, a much better choice is to use the 4th or 
6th order approximation in actual simulations. Still higher order approximations do not 
gain much because they improve the values of the functions only at velocity too close to the 
sound speed. 

Deviations of the stress tensor and of the energy flow are demonstrated in Fig. |21 This 
figure shows that for u ~ 0.075, the gain in exactness of the equilibrium populations greatly 
outweighs the minor deviations in the pressure tensor. On the other hand, Fig. 121 also 
shows that the dominant deviation is in the energy flow. This error is exactly the same for 
all quadratic approximations employed in the standard LBGK simulations, and it can be 
compensated by lowering the viscosity. 

An alternative to the ELBGK model just described, is a straightforward realization of the 
gradient single relaxation time model jl5| . using the H function (fHUJ) derived here. In the 
present context, construction of the corresponding bare collision integral requires only the 
inversion of a 4x4 matrix with populations-dependent entries which can be done analytically. 
This realization does not require any approximation on the equilibrium populations. 

V. DISCUSSION 

In this paper, we derived the H function for the D2Q7 lattice. This makes it possible to 
derive and implement the entropic lattice Boltzmann scheme for the triangular lattice, in 
addition to already established models on square lattices. 

The goal of the entropic schemes is to achieve nonlinear unconditional stability in lattice 
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I III 

10~ 4 10" 3 10~ 2 10~ 1 

~(k) 

FIG. 1: Deviations of the approximate equilibrium population /q from the exact solution (see 
Eq. for /q 9 ). Function e = 10 5 x (/q 9 — /q )/p is plotted for three different values of k. Notice 
that for k = 2, which correspond to the standard second-order polynomial equilibrium used in the 
LBM simulations, the error starts to increase rapidly already at u ~ 0.001. All the quantities are 
given in dimensionless units. 



Boltzmann simulations through creation of valid kinetic models. We recall that the notion 
of the kinetic model of the Boltzmann equation includes the H theorem as one of the 



important properties 



2l| . In the construction of entropic lattice Boltzmann models, 



the second-order polynomial approximations to the generically non-polynomial equilibria of 
the pertinent H functions do appear in the derivations in the same sense as they appeared 
in the lattice gas model, that is, to establish theoretically the hydrodynamic limit to an 
appreciable degree of accuracy in terms of the Mach number. However, these low-order 
polynomial approximations do not show up explicitly in the numerical simulation. Of course, 
this analogy should not be misinterpreted, the entropic lattice Boltzmann equations (0) are 
mesoscopic kinetic equations rather than a lattice gas. As to the numerical efficiency, for the 
already existing ELBGK model, with all the additional burden to solve for the discrete-time 
entropy estimate, the serial processor realization requires only 5 to 10 percent more CPU 
time as compared to the usual polynomial LBGK on the 2DQ9 lattice. 
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10~ 4 10" 3 10~ 2 10~ 1 

FIG. 2: Deviations of equilibrium higher order moments M. Functions e = 10 5 x (AM ELBM )/p 
are shown, where, the deviation from the continuous case is denoted by AM ELBM and are the same 
as the under-braced term in Eq. Q41|) and Eq. (|42j) . All the quantities are given in dimensionless 
units. 

Our final comment concerns the so-called multiple relaxation times models j^. The 
idea behind this approach is as follows: If one uses a second-order polynomial approx- 
imation to the equilibrium, then not all of the linearized collision integrals of the form 
Aj = Y^jLi Aij{fj — ff®) have the same spectral properties, and one can make use of this 
to enhance linear stability by choosing an appropriate matrix Aij. This is done upon con- 
sidering spectra of space-dependent problems (in particular, in a periodic domain 
Although the choice depends on the boundary conditions in the specific spectral problem 
used to determine A^, it might perform better than the standard LBGK also in other flow 
situations. 

To conclude, it is possible to obtain the H function, and a good approximation to the 
correct equilibria for hydrodynamics on the D2Q7 lattice. It was shown, that the quadratic 
polynomial form of the equilibria used in the lattice Boltzmann method is a good approx- 
imation to the correct equilibria for velocity u ~ 0.001. It is possible to obtain a good 
approximation to the equilibria for velocity up to m ~ 0.1 by taking 6th order approxi- 
mation to the correct equilibria. Further, lattice Boltzmann simulations on D2Q7 lattice 
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should avoid using average velocity larger than u ~ 0.075 due to the dominant errors in the 
heat flux. Finally, quadratic polynomials are just a good approximation to the equilibria for 
u ~ 0.001 and should not be confused with the correct equilibria. 
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